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■*->«. I Lattice gas and lattice Boltzmann methods are recently developed numerical 

bJQI schemes for simulating a variety of physical systems. In this paper a new lattice 

r^ \ Boltzmann model for modeling two-dimensional incompressible magnetohydrodynam- 

2 I ics (MHD) is presented. The current model fully utilizes the flexibility of the lattice 

Boltzmann method in comparison with previous lattice gas and lattice Boltzmann 

X" 

Vh ' MHD models, reducing the number of moving directions from 36 in other models 

to 12 only. To increase computational efficiency, a simple single time relaxation 
rule is used for collisions, which directly controls the transport coefficients. The bi- 
directional streaming process of the particle distribution function in this paper is 
similar to the original model [ H. Chen and W. H. Matthaeus, Phys. Rev. Lett., 
58, 1845(1987), S.Chen, H.Chen, D.Martinez and W.H.Matthaeus, Phys. Rev. Lett. 
67,3776 (1991)], but has been greatly simplified, affording simpler implementation 
of boundary conditions and increasing the feasibility of extension into a workable 
three-dimensional model. Analytical expressions for the transport coefficients are 
presented. Also, as example cases, numerical calculation for the Hartmann flow is 
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performed, showing a good agreement between the theoretical prediction and numer- 
ical simulation, and a sheet-pinch simulation is performed and compared with the 
results obtained with a spectral method. 

PACS numbers: 52.30.-q, 52.65.+Z 



1 Introduction 

Lattice gas automata (LGA) methods [1-5], based upon dynamics of cellular automata 
(CA) have attracted considerable attention during the last several years for both 
modeling physical phenomena and simulating linear and nonlinear partial differential 
equations. The lattice gas method is similar to traditional molecular dynamics in 
that a particle representation is employed for microscopic processes such as particle 
collision and streaming, but the dynamics is much simpler. The fundamental idea 
underlying the lattice gas approach is that simple microscopic dynamics may lead to 
macroscopic complexity. 

The first lattice gas automata model was introduced by Frisch, Hasslacher and 
Pomeau (FHP) [0 in a hexagonal lattice for simulating two-dimensional (2D) hy- 
drodynamics. The basic dynamical model comprises particles scattering and moving 
in discretized space and time. The approach of this system is inspired by classical 
statistical mechanics treatments of systems such as the Ising model and simple cellu- 
lar automata models. Although intuitively appealing, the lattice automata method 
for fluids requires an averaging process in order to obtain the macroscopic fluid vari- 
ables and their dynamics, due to the high levels of noise naturally present in the 
discretized particle representation. More recently there has been a trend towards us- 
ing the lattice Boltzmann (LB) scheme instead of the lattice gas automata method. 
Unlike the lattice gas method in which one keeps track of each individual particle, in 
the lattice Boltzmann approach we are interested only in the one-point distribution 
function. While retaining the advantages associated with parallel implementation of 
lattice automata, the LB method is more efficient and accurate computationally, and 
essentially noise-free. 

The history of lattice gas and lattice Boltzmann models for magnetohydrodynam- 



ics (MHD) can be very briefly summarized as follows. The first attempt to model 2D 
MHD with a lattice gas automata scheme was carried out by Montgomery and Doolen 
0, . In their model the basic FHP model is extended to include additional degrees 
of freedom to account for the vector potential. To update the dynamics, some space 
average quantities need to be evaluated. By doing so, the essential feature of locality, 
that characterizes lattice gas systems, is lost. In addition, because of the vector po- 
tential representation of magnetic field the model is intrinsically two dimensional. It 
is noted that the recent lattice Boltzmann MHD model by Succi et al. has a sim- 
ilar limitation. Another MHD lattice gas automata model with pure local operations 
was proposed by H. Chen et al. P, |^. To account for the Lorentz force, the model 
introduced a tensor (i.e., two- indexed) particle representation, and a random walk 
(or, "bi-directional streaming") mechanism. For each one of these particles, there 
are two vectors attached, representing the momentum and magnetic field vectors. 
During the streaming procedure, particles move along one of the two possible vector 
directions with a probability deduced by requiring that MHD behavior is obtained 
macroscopically. Later, S. Chen et al. [|T^ extended the lattice gas automata model 
into a lattice Boltzmann model. The simulation of the two dimensional LGA and 
LB models for problems with free boundary and simple wall boundaries, including 
the two dimensional Hartmann flow and two dimensional magnetic reconnection, has 



achieved reasonable success |10, 11| in test problems 



The random walk MHD LGA model [S, |9|1, and its extension to the LB scheme 



10|, however, have two major problems. First, because of the random streaming, im- 



plementation of wall boundary conditions becomes complicated, requiring increased 



computational memory and computational work |10|. Second and most important 



although both the lattice gas automata model [H, P] and the lattice Boltzmann model 



|10| can be formally extended into three dimensional (3D) space, real 3D implementa- 



tion is impractical due to memory requirements. In order to include a correct Lorentz 
force, for a lattice with A^ moving directions, N x N particle states are needed. In 
3D LGA and LB models, a face-centered-hyper-cubic (FCHC) lattice of A^ = 24 is 
usually employed. For this case, the random walk model needs at least 576 states, re- 
quiring about 1.2 gigabytes for a system of 64^. Thus, to some degree, its actual value 
as a computational tool is diminished because of the requirement of vast amounts of 
memory. 

In the present paper we introduce a new lattice Boltzmann model for MHD that 
requires considerably less memory than previous models, while continuing to offer 
the computational efficiency of the LB approach. The new model is, in essence, 
a reduction of the previous MHD LBE model |]10[, including a smaller number of 



allowed states, while maintaining the symmetries and most of the desirable analytical 
and computational properties of the earlier method. The new model utilizes a "13 
bit" representation on a 2D hexagonal lattice, in contrast to the earlier requirement 
of a 37-bit 2D model. 

The paper is organized as follows: In Section II we will describe the model and 
show how ideal MHD is obtained in the fluid limit. In Section HI the next order terms 
in the Chapman- Enskog expansion are included, adding dissipative effects (viscosity 
and resistivity) to the model. Following this, two sections are devoted to numerical 
tests. Section IV discusses the linear Hartmann flow problem as a first test of the 
model, for several Hartmann numbers H, which parameterizes the solutions. Section 
V describes use of the model to solve numerically for the evolution of the MHD 
sheet-pinch configuration, i.e. the dynamics of a highly sheared planar magnetic 
field. The results obtained with the LB run are compared with those obtained with 
a spectral method run, using the same initial conditions. Discussions of the results 
and of the model are presented in Section VI. Finally, in Appendix A some useful 



tensorial relations for the derivation of the model are presented for completeness; 
and in Appendix B we show how the same principle that is applied to obtain MHD 
behavior on the hexagonal lattice, can be extended to the square lattice. 

2 Description of the Model 

The model described here is inspired by the previous random walk MHD LGA (or, 
CA) model [^ ^ and is motivated by the need for an MHD LB scheme that is com- 
putationally feasible. This requires overcoming the two problems mentioned above. 
The model in this section, for simplicity, concerns two dimensional problems, but its 
principle also applies to three-dimensional models. 

For our two-dimensional system, we use the standard hexagonal grid 0. In the 
vicinity of given lattice point at x, six nearest neighbors are located at positions x-|-ea, 
with Ba = (cos(27r(a — l)/6), sin{2n{a — l)/6)), a = 1, . . . , 6. Instead of using 6x6 
particle states as in previous models, we only consider a subset of them. Each state 
is labeled by the pair of indices {a, a). The positive particle distribution function is 
represented by /^ with a = 1, ... ,6 and a = 1,2, where a is defined relative to e^ in 
the following manner; a = 1 corresponds to the direction a + 1 (mod 6), and a = 2 
to a — 1 (mod 6). 

The evolution of the system consists of a sequence of a streaming stage in which 
the distribution /^ is propagated from each cell to its neighbour cells, followed by 
a collision stage in which the distribution at each cell is redistributed according to 
some conservation laws, as we will see below. The propagation part of the evolution 
for our particular model consists of partitioning the particle distribution into the two 
directions associated with the state (a, a), 

/:(x,T)^(l-p)/;(x + e„T + l)+p/;(x + e.,T+l), (1) 



where T corresponds to the discrete microscopic time and p is a given parameter 
which represents the fraction of the distribution function f^ that propagates along 
the a direction. Notice that this streaming procedure improves in two ways the 
random walk used in p, ^ [1^. First, the motion of the "magnetic" portion of the 
distribution function f^ is always "forward" because the streaming parameter p is 
greater than zero. (For the 36-bit model [TO], the distribution function is represented 



by /^, (a = 1, ... 6, 6 = 1, ... 6) and a fragment of /^ moves in the direction sign{pah)Gb 
while the remainder of the distribution moves in the direction e^. Therefore, there are 
states (a, 6) for which the distribution streams in the direction — e^). In addition, in 
the present model the angle between the two directions e^ and &„ is vr/S for all states 
(a, a). These properties are important for imposing boundary conditions. In addition 
to these twelve states, the model includes a IS*'^ state denoted by /o that represents 
the fraction of the distribution function at a cell that does not advect at all. This 
"stopped" distribution introduces additional freedom in the model that allows to get 
rid of undesirable dependence of the pressure on the velocity that had plagued earlier 
LGA and LB models |2|, [g. 

Associated with each state (a, a) we define the local microscopic velocity v^, which 
is equal to the mean velocity at each cell, and the microscopic magnetic field B^ 

v^ = (l-p)e„+pe,, (2) 

B^ = rea + qe^. (3) 

Although in principle the parameters r, q and p are unrelated, we will see later that a 
connection between r and q is set by the dynamical requirements. Notice that unlike 
the velocity, the microscopic magnetic field does not, on the surface, appear to play 
an active role in the evolution of the system. (However, later we will see that this is 
not the case; see the discussion after Eq. ([T9|).) 



The density p is defined as the summation of all particle distribution functions, 
p(x,r) ^ /o(x,T) + J2Y: f:{^,T). (4) 

a=l cr=l 

In (^) the limits of summation on a and a are given explicitly; subsequently we will 
sometimes suppress them for convenience when no ambiguity is introduced. 

The macroscopic velocity and magnetic field are defined as averages of the micro- 
scopic fields v^ and BJJ, 

pv^E^/:. (5) 

pB = 5:BX- (6) 

The kinetic equation obeyed by f^ can be written by combining the effect of 
streaming, represented by (1) with the coUisional effects, denoted by the symbol f2^, 
arriving at, 

/:(x,T) = (l-p)[/:(x-e„T-l) + fi:(x-e„,T-l)] + 

p[/:(x-e<„T-l) + fi:(x-e.,r-l)]. (7) 

In the present model, for simplicity, we assume that the collision operator f2^ has a 
single time relaxation form [p!0|. 



where r is the relaxation time and /„ '•'^'^^ is the local equilibrium distribution function 
depending on the local particle density, velocity and magnetic field. 

A crucial step in the development of a LB method is the selection of an appropri- 
ate single particle equilibrium distribution function, associated with vanishing of the 
collision operator. This equilibrium distribution function has to be consistent with 



definitions (|^), (^) and (H), and in addition has to give rise to the MHD equations. 
A suitable equihbrium distribution function fulfilling all these conditions is given by 
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/o^'^'^^ = P 
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(9) 



.12 + a C 

where C = 2{p^ — p + 1) . Positivity of the distribution function is guaranteed if v and 
B are sufficiently small. 

To derive the MHD fluid model, we next form the continuum kinetic equation by 
Taylor expanding (^ in the limit of low frequencies and long wavelengths 0. The 
result to lowest order is 



dt 



+ v^v/: = fi^ 



(10) 



Equations for the density, momentum transport and magnetic momentum transport 
can now be found by taking moments of (|T0|) : 

dp 



,^+V.(pv)^0, 

^(pv) + v-n(°) = o. 



;ii) 



d_ 

dt 



(pB) + V-A(°) = 0, 



where the momentum flux tensor 11*^°^ = '}2a,a^l^lfa^^''\ ^-nd the magnetic momen- 
tum flux tensor A*^°) = '}2a,a B^ v^/^^"^^) are deflned in a similar fashion as in the 



36-bit model |^, |T0|- Plugging the expression for f^'^^'^\ Eqs. (|) and @, in the above 
definitions, we may compute the fluxes of density, momentum and magnetic fleld. 

In computing the flux tensors, we need to make use of some of the freedom con- 
tained in our deflnition of the equilibrium distribution (|^) and (^), which depends on 
parameters p, g, r, and a. This flexibility in selecting the parameters will be used 
to eliminate certain unphysical terms that would prevent the appearance of MHD 
equations at leading order, and also to obtain other desirable properties. 

First, we note that the presence of terms that mix the directions e^ and e^j ap- 
pears to be necessary to obtain a correctly structured induction equation. Next, the 
structure of the equilibrium permits appearance of an unphysical pressure-like term 
in the induction equation. This can be eliminated by choice of a relationship between 
r and g, namely, 

r = -,i±^. (12) 

Now we turn to some considerations with regard to selection of the value of the 

parameter p. Two interesting properties of this MHD system should be mentioned. 

First, in the limit case of p = 0, for which the streaming is only along the a direction 

in an (a, a) state, incompressible hydrodynamics is recovered, as expected. Second, 

if the streaming parameter p is changed to 1 — p, the roles played by a and a are 

interchanged. To clarify this point, let us insert the relation between r and q into B^ 

to obtain 

Bl = ^[-{l+p)ea+{2-p)e,i (13) 

2 — p 

and recall that ^r'^ = {1 — p) ea + pea- If p is replaced by 1 — p, then 2 — p and 

1 + p (and thus q and r) also interchange their values. This property holds for 

the distribution function f^ as well; in other words, the same macroscopic MHD 

properties are obtained under this exchange, including sound speed, and transport 
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coefficients. The streaming parameter p is constrained to be between zero and one 
for mass conservation. However, we will see in the next section that choosing p = 
1/2 eliminates spurious terms that would otherwise appear in the momentum and 
induction equations at second order in the spatial expansion of the kinetic equation. 

Turning to the parameter q, we note that, unlike p, q can in principle take any 
desired real value with the exception of zero. For simplicity, we chose q = (2 —p)/\/3. 
For this particular value of g, |vj^| = |B^| and the microscopic velocity and magnetic 
field have the same intensity. Magnetohydrodynamics behavior is obtained for either 
q positive or negative. Notice that changing the sign of q will reverse the direction 
of B. Therefore, this feature is associated to the fact that if the magnetic field is 
reversed everywhere the fluid flow is unchanged. This property was already present 
in the 36-bit MHD CA model |[) and in its LB version 0. 

The equilibrium distribution function for these values of p and q becomes. 



f^(eg, ^ .^ ^_ v"-V+B"-B 

Ja 42^ U + 12 3 [ '^ 

+ 2[(e,-v)(e,-v)-(e,-B)(e,-B)] 

2 B^ 

+ ^[(ea-v)(e.-B)-(e,-B)(e<,-v)] + — 



and 



ft'^ = P 



a 4 n 

- -V 



(14) 



(15) 



.12 + a 3 
Using this form of the equilibrium, after some straightforward algebra explicit forms 

for the flux tensors are obtained as, 

n-°^ = - I 6ij + UkUi[AijM - 6ij6ki] - BkBi[Aijki - 2(5^^(5^] > , (16) 

Kj = P i^ii^kj - 6ikSji)ukBi. (17) 

The ideal MHD equations emerge from this procedure, 

P^r + P(v ■ V)v = -V(P + fiV2) + (B • V)B (18) 

ot 

11 



_ + (vV)B = (B-V)v, (19) 

where the pressure P = pC^. Cg is the sound speed of the system and has a sim- 



ple form related to C and a, Cs = w3C/(12 + a), thus being controllable by the 
parameters p and a. 

Now we make a small digression to clarify some points about the microscopic 
properties of the 13-bit model. For the 36-bit MHD model (LGA or LBE), we recall 
that the macroscopic velocity and magnetic field are defined as follows ||, ||, p^ , 



PV = ^ VLabfab = H [(^ ~ \Vab\)Ga + Pab^b] /, 



ab 



a,b 



pB = Y^ Bab fab = Yl i^abGa + qab^b] fab- 
a,b a,b 

The parameter matrices P, Q, and R (of elements pab, Qab, and Vab, respectively) 
involve, in principle, 108 independent scalars. Arguing that some conditions must be 
placed on these matrices to obtain the right MHD behavior, it was possible to reduce 
the number of independent scalars to only six 0. Although we have not explicitly 
imposed such conditions for the case of the 13-bit model, we want to show that these 
properties are already present in the model as it is. 

For the case of the 36-bit model, because the microphysics should be isotropic, 
we expect that rotating e^ and e?, together by a multiple of 7r/3, Uab and Bab should 
rotate by the same amount. This condition implies that the matrices should be 
circulant. Recall that a tensor H is circulant if for c = 1, ■ ■ ■ , 6, Eat = 'Ea+c,b+c;imod6), 
a, 6 = 1, ■ ■ ■ , 6. In addition, the microscopic physics should be mirror symmetric: if 
Ba and efe are interchanged, then the new values for Uat and Bab, should be the mirror 
images with respect to the line that bisects the angle between e^ and ef,. This implies 
that the matrices are symmetric, Sab = '^ba- These two conditions are trivially obeyed 
by our 13-bit model because the matrices P, Q, and R in the 36-bit model are the 
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counterpart of the scalars p, q, and r. 

There is another constraint to be enforced on the coefficients that is more subtle, 
and is associated with the vector nature of the velocity, and pseudovector nature of 
the magnetic field. There exist microscopic transformations that reverse the direction 
of one of the fields (v or B) everywhere, while leaving the other one unchanged. 
Imposing such a property guarantees, for example that if B ^ — B the evolution 
of the velocity field will remain unchanged, while — B becomes the solution for the 
induction equation. This property was also included in the 36-bit model, by imposing 
some constraints on the parameter matrices, namely, Pab = ~Pab+3, Qab = Qab+s, and 
Tab = —rab+3, where all the sums are modulo{6). It can be easily seen from the 
definitions u„5 = (1 - IPabD^a + Pab^b and B^^ = r„be„ + qabeb for the 36-bit model, 
that by changing every 6 by 6 + 3 Uab is reversed, whereas Bq^ is unchanged. The 
opposite is true if every a is replaced by a + 3. In summary, 

Bafe+3 = ^ab ^ab+S = —'^ab 

Ba+36 = — Bafe Ua+3b = ^ab- 

This property should also be obeyed by the 13-bit model we are dealing with in this 
section. There is no straightforward "translation" from the a + 3 or 6 + 3 operation 
in the 36-bit scheme to our model here, because now for every a there are only two 
6's: a + 1 and a — 1 modulo{6) . Nevertheless, we find that the 13-bit model obeys 
the following relationships, 

Ba+4,2 = BqI Ua+4^2 = —^al 

Ba+1,2 = ~Bal ^a+1,2 = UqI, 

and therefore, such transformations are already embedded in the model as it is. 

As a final comment, we note that this scheme would not exhibit MHD behavior 
in the context of a CA-type lattice gas model. If the idea of the splitting of the 
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distribution in two parts of the present model is "translated" to the CA realm, we 
would, most likely, end up with the 36-bit lattice gas scheme of Refs. p, |[ that 
inspired our model, in the first place. Although this final statement would be hard 
to rigorously prove, we suspect that the present lattice Boltzmann scheme allowed us 
to get rid of all the "degeneracies" (or most of them) present in the 36 bit model. 

3 Transport Coefficients 

In this section we examine in detail the structure of the present model from the 
perspective of a Chapman-Enskog expansion procedure. This renders the long wave- 
length, low frequency behavior of the system, including corrections to the ideal equa- 
tions in the form of dissipative transport effects. We start from the discrete kinetic 
equation (^ 

/:(x,T) = (l-p)[/:(x-e„T-l) + n:(x-e„T-l)] + 
p[/:(x-e.,T-l) + fi:(x-e.,T-l)], 

and expand up to second order in time and space variables, to obtain 

1 do.'' 1 (9^ f 

--[(1 -p)e,e, + pe.e.] : VV(/: + ^l) + ^ - ^^ = 0- (20) 



Now we adopt the following multiple scale expansion [^ |T5[. The time derivative is 

expanded as 

d d r, d ,^,, 

where e is the expansion parameter assumed small, implying that ^2 is a slower time 
scale than ti, and will be associated with diffusion effects. Likewise, the one-particle 
distribution function is expanded, assuming small departures from equilibrium, 
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where /^ *•''•* = fa^^"^- Finally, for the collision operator we write 



r 



(23) 



Replacing all these expansions into pO), we find the following relation to order e, 



dfa 



a(0) 



dh 



+ < ■ v/: 



a{0) 



J a 



a{l) 



T 



(24) 



and 



dh dt2 " ^" " dti'' 

_i A r(i) _ l^Mo) = _1 p(2) 

rdti" 29ti2^" r^" 

to order e^. From ( p^ we can obtain the auxiliary relationship 

1 



(25) 



— (— + v'^■V)r^^)- 






<7(0) 



that can be combined with (|25|) , and after some algebraic manipulations the following 
equation is obtained: 



dt2 ^ 2r' 



dh 



+ K ■ ^)fl 



a(l) 



p(l-p) 



e,-e.)(e,-e.):VV/: 



a(0) 



.1^(2) 



(26) 



Summing equations (^4]) and (p6|) over all velocities, and using that 

'}2a,a fa^^^ = 0, and Y.a,cr fa^"^^ = 0, the following continuity equation up to second 

order is obtained, 



I + V-(pv) = e^^^ E(ea - e.)(e. - e.) : VV/:^. 



(27) 
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The term on the right hand side can be calculated using the tensorial relationships 
(|71|) in Appendix A, and the expression for fa^^'^\ 



|,V,pv)^.^|- 



6 \ dp 3 d{pv^ 



12 + aj dxi C^ dxi 



2p2 _ 2p - 1 
C2 



\, d{pu.u,) ^ djpB') _ ^d{pB,B,] 



dxj 



{21 



dxj dxi 

Therefore, there are second order corrections to the continuity equation. The most 
important term is apparently the one associated with the density diffusion, compared 
to the other terms that are quadratic in the fields u and B. In the limit of hydrody- 
namics, i.e. p = these additional terms vanish. No value of < j9 < 1 (necessary 
for mass conservation) will make these spurious terms vanish. However, notice that 
the r.h.s. vanishes when integrated over the whole domain, and mass conservation is 
restored. 

Similarly, adding moments of equations ( ^^ and ( PB[ ) with respect to v^ and B^, 
the following momentum equation and induction equation are obtained: 

5(pv) P(l-p) 



dt 

dipB) 

dt 



vn 



+ V-A 



$:v^(e,-eJ(e,-e.):VV/: 



a(0) 



p(l-p) 



where 



and 



n 



5:B:(e,-e.)(e,-e.):VV/: 



1 



,7(0) 



(29) 
(30) 



E« /:^°^ + <i-^)/:^'^ 



We can see that there will be several contributions to the transport coefficients coming 
from (|29| ) and (^Of). Let us first examine the contribution coming from the right hand 
side of the equations. For both the momentum and the induction equations only the 
terms linear in the fields in f^^^^ will be different from zero, due to the microscopic 



16 



symmetries already embedded in the model at microscopic level. From the previous 
section we can write 



J a 



^(0) 



+ T;7^[v^v + B^B] + 0(^;^5^). 



12 + a 3C" 
The following tensorial relationships are needed first, 



(31) 



3 3 

Z!(^a)i(ea - e^)j(e„ - e^)k{-v'^)i = ji^ - 3){5ij5ki + 5ik5ji) + -(C + 3) 5u5jk 



5:K),(e,-e.),(e,-eJ,(B:: 



3^3 



{2p - l){6ij6ki + SikSji - 6ii6jk) (32) 



5:(B:),(e„ - e.),(e, - eM^Di = -jiC - 3)((5.,4^ + 

9 

SikSji) + -{C - l)6ii5jk. 



Using the above relationships and after some algebra, we obtain 



p{l-p) 



^v:(e,-e.)(e,-e.):VV/: 



a{0) 



P^^ [2{C - 3)V(V-(pv)) + {C + 3)V^(pv) 
+\/3(2p - 1) [2V(V-(pB)) - V^pB 



(33) 



and 



p{l-p) 



^B:(e,-e.)(e,-e.):VV/: 



a(0) 



V3p{l -p) 



8C 



(2p-l) 2V(V-(pv))-V2(pv) 



+ 



v^ 



-2(C - 3)V(V-(pB)) + 3(C - l)V'(pB) 



The other contribution, that is controllable through r, comes from 



(34) 



n 



(1) 



(1) 






;i-:^)E(B:).(va/:«, 



(35) 
(36) 
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with f^W 
equation, 



-r [d/dt + v^ ■ V] /a *-°'' . We obtain as a contribution to the momentum 

C. 



v-n(i) 



1, 



r 



,C_ 
^4 



+ C,^)V[V-(pv)] + |v^(pv)} 



(37) 



whereas the contribution to the induction equation will be 

i.c r „., 3. 



V-A« = (r - -)^ {-V[V-(pB)] + ^V^(pB; 



(3^ 



We can write the macroscopic equations obtained, including all the above contribu- 
tions to the transport coefficients and in the limit of low Mach number: 

p-^ + p(v ■ V)v = -V(P + pfiV2) + p(B ■ V)B + BV-(pB) 



+ 
+ 



(^-5'|-^'--^' 



V^(pv) 



^-^)(t + ^.') + %7^(^-3) 



V3p(l-p)(2p-l 



V[V-(pv)] 



C 



2V(V-(pB))-V^(pB) 



(39) 



and 



rfVK V 

+ (v ■ V)B = (B ■ V)v + -V-(pB) 
at p 



+ 






(r--)- + ^^^(C-3) 
^2^4 4C ^ ^ 



V^(pB) 

V(V-[pB)] 



v^p(l-p)(2p-l) 



8 



C 



2V(V-(pv))-V2(pv) 



(40) 



where P corresponds to the mechanical pressure. We can readily notice that the 
unphysical terms V(V-(pB)) and V2(pB) in (H), and V(V-(pv)) and V2(pv) in (g) 
can be eliminated in the hydrodynamics limit (p = 0), but more importantly they can 
be removed for p = 1/2, so that the MHD macroscopic behavior can be maintained. 
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The extra terms in the macroscopic equations for p, v and B are also present in 
the 36-bit model described earlier in this work. The origin of these terms lies on the 
bidirectional streaming used in these models. The most seriously offending terms in 
(p9| ) and ( ^OD can be eliminated with the symmetrically appealing choice of splitting 
the distribution in halves. As mentioned in the previous section, the model should 
be symmetric about p = 1/2 and it is noticed that all coefficients appearing in the 
macroscopic equations are indeed invariant under the exchange p —* 1 — p. 

For this choice of J9 = 1/2 (implying C = 3/2) the values for the transport 
coefficients are : 

3r /I 3 \ 1 / 9 \ 

"' = t(4-^I2T^)-4('-^I2T^)- (*2) 

Ai. = -^(3^-2). (44) 

where z/^ and ^h are the bulk viscosity and the bulk resistivity, respectively, and a is a 
free parameter introduced in the previous section that is used to set the sound speed. 
Notice that the ratio u/fi can be arbitrarily chosen by conveniently adjusting the 
parameters r and a. By inspecting these expressions we can realize that the model, 
unlike the hydrodynamics model, displays positive z/ and fi beyond the threshold for 
stability (r = 1/2). Simple stability arguments indicate that the parameter r should 
not be less than 1/2, therefore imposing a lower bound on the transport coefficients. 

4 Hartmann Flow 

We now turn to the application of the lattice Boltzmann model we just described to a 
linear magnetohydrodynamics problem, namely Hartmann flow [W, |T^. This problem 
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represents one of the few MHD configurations that can be analytically solved without 
the need of linearizations (the equations are linear), with the additional assumption 
of constant density and constant transport coefficients. 

The Hartmann configuration consists of a conducting liquid along a uniform chan- 
nel, in steady regime and under the action of a transverse magnetic field. These fiows 
can be used as fiowmeters, by measurement of the potential induced in the fiuids 
as it streams exposed to the external magnetic field ||16|, |18|]. The fiuid, assumed 
incompressible, is constrained to fiow horizontally in a very long, ideally infinite 
channel alongside the x-direction. All relevant quantities, except the pressure, are 
a function of only the transverse coordinate (to the channel) y, v = {vx{y),0,0), 
B = {Bj.{y), Bq, 0). A uniform and time independent pressure gradient is maintained 
along the channel direction to drive the fiuid, so that p = p{x). The walls are located 
at y = —L and y = L. Opposing to the propelling pressure gradient is the viscosity 
of the fiuid and the tension in the magnetic field lines that resist the bending effect 
of the flow. 

For this case, the incompressible MHD equations 

^ + (v . V)v = -V(p + — ) + (B • V)B + uW\, 
— + (vV)B = (B-V)v + /iV2B, 
can be reduced to the following linear system 



i^^ + By^--^ = 0, (45) 



d'^Vx „ dB^ _ Idp 
dy"^ ^ dy pdx 

d^B^ ^ „ dvx 

where we are assuming that the system has reached a steady state; the density p is 
assumed uniform, and By{= Bq from now on) corresponds to the known externally 
applied constant magnetic fleld transverse to the channel. If non-slip (static plane 
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walls) boundary conditions for the velocity, and Bx{—L) = Bx{L) = for the mag- 
netic field are applied, the following analytic solutions can be found for the system 
and (H), 

cosh.{Hy/L) 



"-^^) = i^/!^°'^^^) 



BoP 
BM = ^ 



1 



cosh(if) 



smh{Hy/L) y 



(47) 
(48) 



sinh(if) L 

where the solutions depend on the dimensionless Hartmann number H = BqL/ ^fflU, 
that measures the relative importance of viscous and magnetic forces. In the above 
expressions / = dp/dx, and represents the force driving the fluid down the channel. 
It is easy to prove that in the limit of no external magnetic field Bq (that corresponds 
to if = if all other parameters are maintained constant), the solution for Br,, is 
identically zero, and v^ = —{fL'^/2u)[{y/L)^ — 1]. This is the well-known solution for 
a simple Poiseuille flow. 

The boundary conditions imposed for the magnetic field imply that the walls have 
a finite conductivity (they are neither perfectly conducting nor perfectly insulating). 
In this configuration (with these boundary conditions), the Hartmann flow is operat- 
ing as an electromagnetic flowmeter as we will immediately show |jT7 |. 



The current density j can be obtained from Ohm's law as j = cr(E + v x B), where 
the conductivity a is the inverse of the resistivity /i in our units, and E is the electric 
field. The only surviving component of j is the component normal to the plane of the 
flow, and it can be computed as 

j, = a{E, + v,Bo). (49) 

On the other hand, from Maxwell's induction equation we can get that jz{y) = 
—dB^/dy. This expression for the current density can be immediately integrated to 
obtain the total current across the channel J^. Noting that Bx{y) is an odd function 
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oi y, or simply by using (^) we get 

Jz = - j^ -j^dy = B.,{-L) - B.,{L) = 0. 
Combining with (^) we obtain that 

/L rL 

E^ = -Bo J_ v^{y)dy, 

from which it follows that 

E, = -BoVM, (50) 

where vm represents the mean velocity of the flow across the channel. Consequently, 
vm can be computed by measuring E^. {Bq is assumed to be known since it is 
externally applied.) 

Now we turn to the numerical simulation of Hartmann flow, making use of the 
lattice Boltzmann model with 12 moving states described in previous sections. The 
simulation domain was for all cases a lattice of only 4 cells in the x direction of the 
problem x 60 cells in the transverse direction y. The x direction coincides with the 
direction ei of the hexagons. The system is initialized by setting the distribution 
functions f^ and /o to their equilibrium values given by a uniform density, and a 
transverse magnetic field By ^ 0. B^, Vx, and Vy are initially zero. 

The system is evolved by using the standard sequence of collisions and streaming 
processes, with the addition of an intermediate step that acts to generate an effective 
pressure gradient. To achieve this, the distribution of moving particles at each cell 
/^ is redistributed so that the total density and the magnetic field at the cell are 
unaltered, and the velocity receives a "kick" in the direction of the flow. Care must be 
taken that the redistribution of mass density along the different states does not push 
the distribution function too far out of its equilibrium value, nor makes it negative. 
That is, for the problem under consideration, a forcing function must be constructed 
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that tries to increase pv^ while not changing pvy, pB, or p. To exphcate the dynamics 
of this procedure, let us recall the definitions of the macroscopic quantities, 

p = /o+E/: 

pv = ^[(l-p)e,+pe<,]/; 
pB = Y.[re, + qe^]f:. 

Therefore, if at a certain microscopic time T the distribution function at a specific 
cell is /^ before this "forcing" scheme, and f^' = fa ~^ F*Caa after the distribution 
is kicked, we wish to find a quantity F*Caa with the constraints 

P = T.fa' = T.f: + ^P (51) 

aa aa 

E /:'[(! -P)ea + pe.] = E/:[(l-p)ea+pe.] + p(Av) (52) 

aa aa 

Efa'iirea + qe^] = J2 f^HrSa + qe^] + A{pB), (53) 

0X7 CUT 

where A represents a "small" change. As stated above, we are seeking A(p) = 0, 
A(pB) = 0. Using that f^' = fa + F*Caa and adding over a, Eqs. 51-53 become 

A(p) = = E(C.^ + ^a)=0, (54) 

a 

A(pv) = F*{l-p)Y.{Cl + Cl)ea + F*pY.{Cl_, + Cl^{)e^, (55) 

a a 

and 

A(pB) = = r E(C] + Cl)e, + 9E(Ca-i + ^a+i)e.. (56) 

a a 

Combining (^4|) and (|56|) we obtain 

A(pv) = {l-p-p'-) Y.{Cl + Cl)ea. (57) 

y a 

Choosing A{pVx:) = 1 and A{pVy) = the following set of equations can be found 

T.icl + CD = 

a 
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C2 + C2 + c^ + C^ — C^ — Cr, — Cq — Cq — 



(1-P-P-) 


Cl + Cl + 


Cl , C7| 

2 2 


2 


/^2 /^l /^2 /^l /^ 

•-"a r<i r<2 '-'5 '-'5 , "-"e , '-'f 
2 ^^ ^^ 2 2 + 2+2 


for which a particular solution is 




Cl = 


d = -d = -Cl = 1 




c^ - 
c^ - 


^'~ 2-p 
ri 3 - 2p 



1, (58) 



(59) 
(60) 

(61) 

where we used that r/q = —{l+p)/{2—p). By appropriately choosing F* we 
are able to control the strength of the forcing scheme. We note that the pressure 
gradient produced in this way is uniform, and that it follows the spirit of the CA 
interpretation of cell populations as "particles." When this process is put into action, 
the total momentum in the direction that is being forced is seen to increase, until the 
driving force is balanced by the action of the viscosity of the fluid and the reluctance 
of the magnetic field lines to be bent, and the system reaches a stationary regime. 

The boundary conditions were implemented by setting Vx, Vy and B^ equal to 
zero in the first and last layers (i.e., y = and y = L), combined with a periodic 
streaming. An alternative way to achieve the boundary conditions would be to cancel 
the periodic streaming of the populations, and to make the "particles" bounce off the 
walls reversing the velocity at those cells while keeping unchanged the magnetic field. 

Simulations of the Hartmann flow were carried out for Hartmann number H = 
(zero magnetic field case), 1,3,5,8, and 13. The results are presented in Fig. 1 and 
Fig. 2. Figure 1 displays all the velocity profiles Vx{y) plotted versus y, for all six 
values of H, and compared with the analytical solution (^Tf). A very good agreement 
between the latter and the computational results is obtained for this range of values 
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of H. 

The flattening of the velocity profile can be understood in several ways [0, |1^, ^ 
From the linearized MHD equations (^) and (|4^), it can be seen that the transverse 
magnetic field tries to eliminate vorticity. If the magnetic forces dominate (i.e., for 
large Hartmann number if), then the velocity profile tends to flatten. ||T9| However, 
the velocity profile cannot be flat all across the channel, because the velocity at the 
walls vanishes, so there must be a region within a certain distance from the walls on 
which the gradients of the velocity are very large, i.e., the region where the vorticity 
(produced by the boundaries) is confined to. 

Alternatively, the flattening of the velocity profile could be understood as follows. 
Combining (^Q]) and (|50|) we obtain 

jz = (tBo{v^ -vm)- 

The magnetic force on the fluid is given by the Lorentz force F^, = J x B, that in 
our case reduces to (-Fl)^ = —(TBl{vx{y) — vm), and therefore the flow tends to slow 
down where Vx > vm and tries to speed up where Vx < vm, producing a flattening of 
the velocity profile. 

The high degree of agreement seen for the velocity profile across the channel is 
also observed for Figure 2, that shows the same comparison for Bx{y), for if = 1,3 
and 13 only. For this figure, only results for three different Hartmann numbers were 
included for the sake of clarity, since the effect of H on Bx, is not as marked as it 
is for Vx- We note that the fit for the cases not shown is as good as those displayed 
in the figure. There is a point we would like to stress about the material presented 
in this figures. These are not "fits" in the usual sense. The Hartmann number H is 
here constructed from u and fi which come from the Chapman-Enskog theory. The 
analytic solutions (|47|) and (|48|) depend only upon /, fi, u, Bq, L, and p. The forcing 
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/ is fixed; Bq is fixed initially as are p and the simulation size L. With this subtlety 
in mind, the solutions have no free parameters to adjust. 

Some caution must be exerted when calculating the Hartmann number corre- 
sponding to the simulation, since the exponential character of the solutions would 
make them very sensitive to a small departure from the actual value of H. In partic- 
ular, the width of the channel 2L should be evaluated as (A^^, — l)-\/3/2, taking into 
account the x — y ratio for the hexagons, and the position of the boundary in the 
simulation is at y = 1 instead oi y = 0. 

The Hartmann number H = BqL/ ^/JIi' was varied by changing the strength of 
the magnetic field Bq. The size (60 cells) was kept constant for convenience in the 
manipulation of data. The kinematic viscosity v and resistivity /i can be obtained 
as a function of the relaxation time r. From Section III we recover the expressions 
V = 3r/16 and /i = 9r/16 — 1/4. Only r = 1 was used for all the simulations in this 
section, therefore the system is forced to equilibrium in each iteration. The "forcing" 
coefficient was set to F* = 2 x 10~^, thus being sufficiently small (comparing F*C^ 
with the mean value of the density, po = 3.9) to ensure that the distribution function 
f^ will be only slightly departed from equilibrium during the "forcing" step. 

The limitations for the range of H that the model will be able to accurately 
reproduce are as follows. 1) for the upper bound, we can see from (|^) and (^) 
that we can estimate the width of the boundary layers, in which the Lorentz force 
is comparable to the viscous forces, as 5 ~ y/17JI/Bo, from which we gather that 
6/L ~ H^^. When H ^ 1, the thickness 6 becomes very small (the region of 
nonzero velocity gradients is confined to a very narrow layer away from the walls), 
and chances are that we need to increase the width of the domain to resolve 6. For our 
simulations, for which 60 cells across the channel were used for all cases, we observed 
that for H > 30 the boundary layer thickness 6 is of the order of one cell. Therefore, 
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if simulations with higher values for H are desired, L should be increased. We note 
that this is not a limitation of the lattice Boltzmann simulation scheme, but rather 
a resolution constraint: the computation of the analytical solution presents the same 
flaws. 2) for the lower bound, the limitation is given by the roundoff error of the 
machine since for /J <^ 1 we are forced to use very weak external fields Bq. 

The one dimensional nature of the problem is highlighted by the fact that the 
exact behavior of the solutions was reproduced for a domain that was 4 cells long in 
the fluid direction x. As a matter of fact, we observed that the same results can be 
obtained with a length of only one cell, making the domain truly one-dimensional, 
and supporting the hypothesis made for the derivation of the lattice Boltzmann ap- 
proach, that the population of the cells correspond to ensemble averages of the discrete 
populations used for the Cellular Automaton approach. 

5 2D Magnetic Reconnection 

In the previous section, we argue in favor of our lattice Boltzmann MHD model, 
by comparing its solutions for the Hartmann flow problem, with the analytically 
obtained solutions. The results are encouraging to the extent that is very difficult to 
observe with the naked eye in Figures 1 and 2 departures of the LBE solution from 
the theoretical ones. This was the case for a wide range of values of the Hartmann 
number, the only parameter in the problem. Optimistic as we might be, we recognize 
that Hartmann flow is essentially a one- dimensional and linear problem. Thus, it is 
our intention in this section, to test the validity of the 13-bit LBE model for a situation 
that is both two-dimensional and nonlinear. The configuration we chose is the 2D 
MHD "sheet pinch" . Before presenting results for the present 13-bit model, we recall 
that the reconnection configuration had been chosen previously as a test problem 
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for 36-bit MHD model of Ref. |[T0[- The authors found that the LBE solutions 
were qualitatively correct, and they showed features of the evolution associated with 
nonlinear effects of the sheet pinch dynamics. However, no claim of comparison with 
the "real" solutions, or solutions from other numerical methods was made. 

The "sheet pinch" is a magnetohydrodynamic configuration characterized by an 
inhomogeneous magnetic field that changes markedly in a very narrow region, thus 
producing very strong currents. This arrangement can be encountered in a variety of 
important physical phenomena as solar flares, and the earth's magnetic field, reversing 
its sign embedded in the solar wind. It is not the goal of this section to discuss in 
detail the physical processes in a magnetofluid undergoing a reconnection process, a 



vast literature is available. Central are the theoretical efforts of Dungey [20|, Parker 



pT| , and Sweet ||22[ . For a discussion of the reconnection process, including the role 



played by the fluctuations from the point of view of turbulence, see Matthaeus and 



Montgomery [Q, and Matthaeus and Lamkin [p^ . 

Our approach to the "sheet pinch" problem will be similar to the procedure fol- 
lowed in another recent test of the LBE method ||25|, in which LBE and spectral 
method solutions for a 2D hydrodynamic shear layer were compared in detail. Here, 
we will briefly describe the reconnection runs from a technical point of view, and 
then we will move on to present and contrast the results obtained with both LBE 
and spectral methods. 

5.1 The Sheet Pinch Simulations 

The idealized sheet pinch consists of a uniform magnetic field reversing sign in a very 
thin zone, much in the same way as the velocity swaps its direction in the idealized 
nonlinear Kelvin- Helmholtz instability. This configuration gives rise to a current sheet 
because j = (V x B)^, with j the current density in the z direction. Therefore, the 
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initialization of the sheet pinch was done, in a 27r x 27r simulation box, with a spectral 
truncated representation of delta functions located at y = 7r/2 and y = 37r/2 (with 
opposite signs) for the current, including wavevectors k = 1 through 31. 

An uninteresting evolution follows unless some non-sheet pinch modes are excited, 
to initiate the nonlinear couplings. The current density and the vorticity Fourier 
modes with 1 < A; < 15 where excited with random phases and with spectra of k~^ 
for high k for both, kinetic and magnetic energy. The "noise" spectrum was peaked 
at A; = 3 for both kinetic and magnetic energy, and added about 1% of the energy 
already present in the ideal sheet pinch. 

For the spectral run the z components of the vorticity u = (0, 0, u) and the vector 
potential a = (0, 0, a), are evolved according to the equations 

-^ + V ■ Vcj = B ■ Vj + z/V^o^ (62) 

da n 

_ + vVa = /iV^a, (63) 

at 

where fi is the resistivity and u is the kinematic viscosity. The fields v and B are in 
the plane x, y depending solely on those coordinates. The scalar functions uj and a 
are related to v and B by c<jz = V x v, and B = V x az = Va x z, where z is the unit 
vector in the direction normal to the x,y plane. The current density is V^a = — j, 
and the vorticity is related to the stream function by V^ip = —uj. The spectral run 
is of the Fourier-Galerkin type and has a resolution of 128 x 128. 

For the LBE run we have to specify the initial density, velocity field, and magnetic 
field. The velocity field is obtained from the relation v = Vip x z, the stream function 
ip being obtained from the solution to V'^ip = —uj, with uj being the initial vorticity 
used in the spectral run. Similarly, to initialize B, we numerically evaluate B = Vaxz 
in Fourier space, a being the initial vector potential used for the spectral run. The 
initial density is set to a uniform value p = Po = 3.9. The initial fields are used 
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to evaluate the equilibrium distribution function /^ for this model, given by Eqs. 
(H) and d), for these initially specified fields. From then on, the LBE sequence of 
streaming and collisions is performed to evolve f^. 

A short digression is needed at this point to justify the choice of z/ and fi for the 
spectral run (the inverse of the Reynolds number R, and the magnetic Reynolds num- 
ber Rm, respectively, in the set of units we are using), and the relaxation parameter r 
and the simulation size, for the LBE simulation. Since it is our goal to test the MHD 
LBE model in a situation that involves strong nonlinear interactions, it is desirable 
to perform a simulation with the Reynolds numbers R and R^ as large as possible. 
The limitation is imposed by the LBE scheme. In Section III, explicit expression for 
the resistivity and viscosity of the LBE model were found: 

3 

i^LBE = —r (64) 

lo 

9 1 

For stability reasons, the parameter r is constrained to be r > 1/2. On the other 
hand, for attaining small i'lbe and fiLBE, and thus high R and R^, the relaxation 
parameter r should be as small as possible. Choosing r slightly larger than 1/2, to 
ensure stability, the values of R and Rm will be dictated essentially, by the simulation 
size. To strike a compromise between a computationally reasonable simulation size, 
and the degree of nonlinear activity, we use a 512 x 512 resolution domain, for the LBE 
simulation. The characteristic speed in our problem is the Alfven speed, coinciding 



with J {B"^) in our units. Consequently, the magnetic Reynolds number is given by 

BL 512 1 

Rm = = 0.1x^— X-, 66 

/i 2-11 fi 



where B = J {B'^) = 0.1 initially. Similarly, for the mechanical Reynolds number we 
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obtain, 



fi=^.0.1x5Hxi. (67) 

V 271 V 



For r = 0.5001 we get R = 86.9 and Rm = 260.3. The reciprocal of these numbers 
are those used for z/ and /i, respectively, for the spectral run. The factors of 2tt are 
necessary because the length L used to define R and Rm for the spectral simulation is 
the physical length of the simulation box, divided by 27r, i.e., L = 27r/27r = 1. Thus, 
for the LBE simulation simulation we need to use 512 (the physical size) divided by 
27r. Note that for the LBE run, because we are using the hexagonal lattice, the typical 
length used for the above evaluations, is somewhat ambiguous due to the V^/2 ratio 
between the y and x directions. Moreover, the system represented by both methods 
are physically different. 

Last, before turning to the comparison of the results obtained for the two runs, 
and to be able to relate physical processes observed in both simulations, we need to 
find a relationship between the spectral and LBE characteristic times. Let us write, 

Tlbe _ Llbe/Lsp _ 512 1 

~T^ ~ Ulbe/Usp ~ ^OT' ^ ' 

where Ulbe = 0.1 is the characteristic speed of the LBE model, given by the Alfven 
speed, and similarly Usp = 1 for the spectral method run, that tells us that Tlbe = 
814.87 Tsp. 

The runs were carried out up to about ten spectral-method characteristic times. 
Periodic boundary conditions were imposed for the simulations. 

5.2 Comparison and Discussion of the Sheet Pinch Runs 

The evolution of the sheet pinch dynamics, being an MHD system, is much more 
complex than its hydrodynamics counterpart. There are more dynamical variables 
evolving coupled to each other, and there is a larger parameter space. For example, 
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the relative values of R and Rm may influence the dynamics. Nevertheless, there are 
some global features thought to be similar for most decaying, 2D, incompressible MHD 
flows in periodic geometry. For example, there are three rugged invariants |^, that 
is constants of the nondissipative evolution that survive the truncation in k space 
(Galerkin approximation). Those are ||2^, the total energy {E) = X]k[^^(k)/A;^ + 
a'^{'k)k'^], the cross helicity H^ = {v -B) = J^k^O^) c^lk), and the mean square vector 
potential A = {a?) = J2k(^'^(X)y where (■ ■ ■) denotes a volume average. 

In Fig. 3 we present time histories of bulk quantities that characterize the turbu- 
lence for both the spectral method and the LBE run. The continuous line corresponds 
to the spectral method simulation. Panels a) and b) of Fig. 3 display the evolution 
of A and E, respectively. Both quantities, which are ideally conserved, decay mono- 
tonically due to the presence of nonzero dissipation coefficients /i and z/. The slow 
decay of A suggests that there is a dynamic redistribution of this quantity in favor of 
larger scales. In Fig. 3c) we present the evolution of the kinetic energy. The kinetic 
energy Ek, initially 1% of the total energy, decreases throughout the run, displaying 
some small bursts of activity that would be more pronounced for larger R and Rm 

In panels d) and e) of Fig. 3 , we show the evolution of the enstrophy Q and 
the mean square current J = (j^). These two quantities highlight the activity in 
small spatial scales, and for these runs they rapidly decay due to the relatively high 
viscosity and ohmic dissipations. Again, both quantities present a more "bursty" 



shape for more turbulent systems |24 



We now turn our attention to comparison of contour plots. In Fig. 4 we exhibit 
plots of constant magnetic field lines (constant a), for times approximately equal to 
one, seven, and ten. Emergence and subsequently growth of magnetic islands can be 
seen. For the same times, in Fig. 5 , contours of constant vorticity are displayed. 
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Although the distribution of vorticity consists of small, nonlocalized structures, it 
can be seen that this quantity rapidly organizes in the region of the X-points to form 
quadrupole-like structures. The vorticity plots, combined with the stream function 
contour plots, shown in Fig. 6 , provide a consistent picture of part of the activity 
of the magnetofluid in the reconnection zone. Jets of fluid are seen to come into the 
"hot" area (higher speeds are represented by denser ip lines) from the strong-field sides 
of the X (up and down in our plots), and out through weak- field corners (sides of the 
X). These jets are responsible for the bursts that suggest themselves in the evolution 
of the kinetic energy. This is essentially a pressure-driven effect due to the steep 
gradients of the magnetic field present near the neutral sheet; the fiuid finds itself 
pushed by the magnetic pressure into the neutral zone, and, being incompressible, 
it has no "choice" but to turn into the weak-field region producing the four vortices 
seen in the plots, at the X points. 

The last set of plots from the spectral and LBE runs (Fig. 7), shows contours of 
constant current density. Initially the current is concentrated, more or less uniformly, 
along both current sheets. As the system evolves, and the magnetic field lines start 
to reconnect, we see a tendency for filaments of current density to form. The regions 
with current filaments will participate in a significant part of the energy dissipation 
due to finite resistivity. 

We conclude from the examination of these plots, that this hexagonal, two- 
dimensional MHD LBE model, is capturing the basic mechanisms of MHD dynamics. 
From the contour plots, we see that once a structure of one of the quantities is identi- 
fied in one of the plots corresponding to the spectral run, a very similar structure can 
be also seen in the corresponding LBE plot. Similarly, the LBE tracks the evolution 
of relevant bulk quantities very closely. We would like to point out that a perfect 
agreement is not, in fact expected. Additional refinements could be introduced in 
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the LBE simulation to still further reduce the gap between the results from both 
simulations. For example, it was mentioned above that the spectral simulation box, 
and the LBE domain are physically different systems, due to the effects coming from 
the hexagonal lattice. This will certainly affect the magnitude of the mechanical and 
magnetic Reynolds numbers. There are, at least, two alternative ways of getting 
around this difficulty. One of them would be to make use of the similar MHD model 
derived based on the square lattice, which is described in Appendix B. (The earlier 
comparison of the LBE and spectral method hydrodynamic shear layer dynamics ^E 



employed an LBE on a square lattice.) Although the hexagonal LBE requires less 
memory, this is not a decisive advantage. Another possibility would be to abandon 
the use of the same number of cells in both dimensions {N^ = Ny), and to choose, 
for example Ny = Nx\/3/2. This choice would have introduced complications in di- 
agnostics currently based upon Fast Fourier Transforms. Instead, considerable larger 
amounts of data would be required to be kept for later analysis. 

An examination of the divergence of the magnetic field is mandatory, to make 
sure that monopoles are not being created by the model, thus casting doubt on the 
results. To this end, we decompose in Fourier space the magnetic field B(k) in its 
longitudinal component Bl, and its transverse component B±. A similar examination 
was carried out for the velocity field in the shear layer LBE study [§^], as a way to 
quantitatively measure the compressibility of the flow. We calculate B^ and B± as 

Bl ^ ^^^2imi. po) 

This ratio is a good measure of the amount of "monopolar" (longitudinal) activity as 
compared with the transverse component, containing most of the energy. In Fig. 8 we 
show the evolution of the ratio Bl/B±. We readily notice that the overall tendency 
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of this quantity is to decrease, and that by the end of the run the "amount" of Bi is 
about one part in one thousand. 

The nonzero initial value of B^ is attributed to the non-square nature of our 
LBE simulation box. Exactly the same field that produces V-B = for the square 
spectral simulation run, produces nonzero divergence on the hexagonal lattice. What 
is encouraging about this picture is that the LBE dynamics seems to possess self- 
adjusting mechanisms that reduce the amount of monopoles, much in the same way 
as Chen et al. |^ discuss in the context of their 36-bit MHD CA model. 

Finally, we turn to a brief discussion of the efficiency of the method. At the end 
of the previous section, we noted that although the transport coefficients are directly 
controllable via the relaxation parameter r, the threshold of stability with respect to 
T, is higher than the value of r needed to make v and /i zero. This technical problem 
limits the Reynolds numbers that can be attained in the MHD LBE for fixed grid 
size. In the case of the hydrodynamic shear layer LBE |^, the simulation domain 
was chosen with the objective of resolving the spatial structure of the turbulent ac- 
tivity, much in the same way as it is done for a spectral method simulation. Thus, at 
the same Reynolds number, an LBE and spectral method hydrodynamics simulation 
can have about the same grid size. For the MHD LBE, the size was also determined 
by the requirement of matching the Reynolds numbers with the spectral code. This 
512 X 512 LBE simulation was run on a CRAY-YMP computer, in the San Diego 
Supercomputer Center, and needed about 12 minutes per characteristic time (about 
800 LBE microscopic times). The spectral run, 16 times smaller (128 x 128), took 
about ten times less CPU time. Nevertheless, it should be noted that the numer- 
ical efficiency of LBE-type computations is greatly enhanced in massively parallel 
computers, due to its local dynamics. 
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6 Discussion and Conclusions 

In this paper we have introduced a model for simulation of 2D MHD with the lattice 
Boltzmann equation technique. The idea of propagating the distribution at a given 
state into two directions associated with the velocity and magnetic fields, had been 
previously used for obtaining 2D MHD using CA dynamics [||, |^ , and later extended 
to LBE [l^. In the present scheme, by utilization of the same idea combined with the 



flexibility of the LBE scheme a significantly more efficient and simpler method is ob- 
tained for simulation of 2D MHD. The improvement is two-fold; first, the number of 
discrete velocities is reduced in our model from 37 to 13, only. Second, the algorithm 
for the evolution of the present model is simplified by requiring a "forward" stream- 
ing, as explained in the text. These models possess the same microscopic symmetries 
necessary for guaranteeing the correct long wavelength, low frequency behavior. The 
theory for the model was presented, including the Chapman-Enskog expansion pro- 
cedure to obtain second-order effects. In passing, we note that the simplicity of the 
model is apparent when evaluating all the second-order contributions displayed in 
Section HI. 

Evidence of correct MHD behavior was introduced in Sections IV and V, where 
the model is applied to reproduce a linear problem (steady Hartmann fiow), and a 
nonlinear problem (evolution of the 2D sheet-pinch). For the former, the performance 
of our model is extremely good, for a reasonably wide range of the Hartmann number, 
that parameterizes the problem. The second numerical test of the model behaves 
reasonably well. Nevertheless, it should be mentioned that this more stringent test 
exposes what might be the most serious deficiency of the model, namely the inability 
to achieve relatively low transport coefficients. This disadvantage hurts the potential 
use of this model for highly turbulent simulations. Clearly further investigation into 
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this matter is required. At the moment we ignore the fundamental reason of this 
crossover (in terms of the relaxation time r) between the stability threshold (r > 1/2), 
and the region of low values of viscosity and resistivity that occurs for r < 1/2 (see 
Eqs. (|4l| ) and (|^)), and whether this effect is induced by the special streaming used 
for the model, or the deficiency could be cured by choosing a more flexible collision 
operator. The LBE model for hydrodynamics (with single-time-relaxation collision 
operator) does not share this inconvenient feature since r > 1/2 is both a condition 
for numerical stability and for positivity of the viscosity. Another two desirable 
properties that an improved MHD LBE model could have are as follows; first, it 
would be convenient to have independent control on the viscosity and resistivity (in 
the present model, the choice of the relaxation time r determines both /i and u). 
Second, and more important, in the present model, and in other MHD models we 
made reference to in this report, the divergenceless property of the magnetic field is 
not imposed. It is seen, however (see Ref 0) that the model possesses some self- 
adjusting mechanisms that diffuse away the solenoidal component of the magnetic 
field, as displayed in Section V, for the sheet-pinch simulation. 

In spite of these improvable aspects, we believe that the 13-speed model is cap- 
turing the essential features of the equations for MHD with a minimum number of 
degrees of freedom, and that it is a valuable tool for non-turbulent regimes. More- 
over, although the model was formulated explicitly for 2D in the present paper, its 
extension to 3D with a reasonably low number of degrees of freedom does not seem 
to pose any serious difficulty. 
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Appendix A: Useful Tensorial Relationships 

The relevant tensorial relationships used extensively for the derivation of the fluxes 
of density, momentum, and magnetic field in Section II are: 

3 

a,cr a,a ^ 

3 
3 9 

a,(j ^ ^ 

^{ea)i{,Ga)j{ea)k = Xl(e^)i(e^)j(e^)fe = (71) 

I](ea)i(ea)j = ^(e^)i(e^)j = 6% 

a,a a,(T 

IZ(e«)i = J2(^^)j = 0' 

a,cr a.cr 

where Sij is the Kronecker delta, and Aiji^i = Sij6ki + SikSji + SuSjk- 

Appendix B: Simulating 2D MHD on the "Square" 
Lattice 

Moving from modeling hydrodynamics to modeling MHD requires the inclusion of 
new force terms, and more importantly a new whole equation to follow the evolution 
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of the magnetic field. This equation involves terms like (v • V)B — (B • V)v, and the 
technique that has worked in the 36-bit and in the 13-bit models, has been splitting the 
distribution in different directions during the streaming part of the evolution. This 
parting of the distribution produces a mixture of directions that allows to include 
terms mixing v and B and with different signs. Once this fact is recognized, it is not 
hard to generalize the same idea to the square lattice. 

In this section we would like to simply display a method for modeling 2D MHD on 
the square lattice. This would have a slight advantage and a slight disadvantage over 
the hexagonal lattice (13-bit model). The latter is that the "square" model requires 
more memory than the hexagonal model: we need to keep track of 13 real numbers 
per cell, for the hexagons, whereas for the squares the requirements increase to 17 
real numbers per cell. On the other hand, the square lattice provides a "better" 
simulation domain, in a geometrical sense, for some systems, unlike the hexagonal 
case, for which the basic cell has different physical lengths in the two directions. 

We now turn to a brief description of this alternative approach. If one of the 
square cells is located at x, its nearest neighbors are located at the face-centers 
X + c^, for a = 1,2,3,4, with c^ = (cos (a — l)7r/2, sin (a — l)7r/2), and the ver- 
tices of the square centered about x, i.e., x + c^J , for a = 1,2,3,4, with c^^ = 
V2(cos (a - l/2)7r/2, sin (a - l/2)7r/2). 

There will be three distribution functions: one that streams in the lattice indicated 
by the superscript /, a second one that is moved in the lattice //, and a "stopped" 
distribution. Therefore, the streaming part of the evolution can be represented by 
the expression, 

/af(x,T)^i/iI(x + cf,T+l) + i/j^(x + cf,T + l), (72) 

where K = I or II, a = 1, 2, 3, 4, 6 = a -|- 1 or a — 1 (modulo 4), and T is the discrete 
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lattice time. 

The macroscopic quantities we are interested in following, the density, velocity, 
and magnetic field, are defined below. 



P = fo+Efa, 
a,b,K 

^ a,b,K 

pB = E?i^(-cf + cf)/,t 

a,b,K 



(73) 

(74) 
(75) 



where gi = 1/2, and g2 = 1 and /o represents the stopped distribution. 

During the collisional part of the evolution is when the three types of distributions 
"see" each other. We use the single time relaxation approximation with parameter r, 
so that the discrete kinetic equation obeyed by the system is, 

1 



/j^(x,T) = -[/„f(x-cf,T-l) + f]„,(x-cf,T-l)] + 
^[/„^(x-cf,T-l) + (],,(x-cf,T-l)], 



(76) 



where 



n, 



ah 



■(/< 



K 

ab 



fab^''^)/r. 



(77) 



The procedure to get the macroscopic MHD equations is familiar to us at this point. 
The only thing that is left to complete the definition of the model is to specify the 
distribution functions. 



J ab 



/; 





ab 



pd 



K 



1 + 



1 



24rf, 



(cf + cf)-v + ^(cf + cf) 



32^2 



-K 



+ 



1 



24rfo 



B + 



32^2 



fc5 ■ v)2 



1 



24^2 



c„--v)(cf ■B)-^(c-.B)(cf -v) 
^(cf-v)(cf.v)-^(cf.B)(cf.B) 



pdo 



2dn 



(78) 
(79) 
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where 



do = 1- 40(i2 di = 4(^2 

«! = 1 "2 = 2' 

K = I or II, and < ^2 < 0.025 for positivity of lio- 

For the sake of completeness, we document the tensorial identities relevant to the 
derivation of the model 

Ecfcf = Ecfcf = 2A^I, 

ab ab 

Ecfcf = 0, 

ab 



E(cf ).(cf ),(cf )fc(cf ), = 2ZkA,,,i + 2YKS.,ki, (80) 

ab 

E(cf)^(cf).(cf)^(cf)^ = o, 

ab 

E(O.(cf).(cf),(cf), = 0, 

ab 

E(cf ).(cf ),(cf )fc(cf )i = Al6,,6ki - 2ZKA,,ki - 2YK5,,ku 

ab 

where Aijki = Sij6ki + 6ikSji + SaSjk, 6ijki = 1 only Hi = j = k = I, otherwise is 0; 
Ai = 2, A2 = 4, Zi = 0, Z2 = 4, Yi = 2, and Y2 = —8. / represents the identity 
matrix. 

The viscosity and resistivity can be calculated using the Chapman-Enskog expan- 
sion, in a similar fashion as it was done in Section III, 

r + 1 , , 

^ = — (81) 

A = ^. (82) 

The transport coefficients have been numerically measured for the case of decaying 
shear flows, and the values found were in agreement better than 0.1% with the pre- 



dictions of (|81|) and (H). 
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Figure Captions 

Fig. 1 Profiles of v^ vs. y/ L for Hartmann number iJ = 0, 1, 3, 5, 8, and 13, 
from top to bottom. The H = (unmagnetized) case corresponds to the 
Poiseuille flow. The fines indicate analytical results, and tfie symbols are 
tfie solutions provided by tfie LBE scfieme. 

Fig. 2 Profile of B^ vs. y/L for H = 3 {+ symbol), H = 1 {o symbol), and 
H = 13 {A symbol). Tfie lines indicate analytical solutions, and tfie symbols 
represent tfie LBE solutions. Simulations witfi H = 5 and 8 were carried 
out witfi similar results, and tfiey are not presented fiere for clarity. 

Fig. 3 Time fiistories of bulk quantities for tfie LBE run and tfie spectral run. 
Evolution of a) tfie mean square vector potential (A); b) tfie total (magnetic 
plus kinetic) energy E; c) kinetic energy (Ek); d) mean square vorticity (tfie 
enstropfiy Q), and e) tfie mean square current J. Tfie LBE run is indicated 
by tfie solid line and tfie spectral run by tfie dasfied line. 

Fig. 4 Contours of constant magnetic field (constant a), for tfie spectral (SP) 
and LBE runs, at times approximately equal to one, seven, and ten. Growtfi 
of magnetic islands can be observed for botfi runs. 

Fig. 5 Lines of constant vorticity at times approximately equal to one, seven, 
and ten, for botfi runs. Quadrupole-like structures can be noticed in tfie 
region of tfie X-points. 

Fig. 6 Contours of constant stream function (ip) at times approximately equals 
to one, seven, and ten. Similar features can be observed for botfi spectral 
(SP) and LBE runs. 
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Fig. 7 Lines of constant current density at times approximately one, seven, and 
ten, for the LBE run and the spectral run. Diffusion from the sheets area 
can be observed, as well as filamentation in the X-point regions. 

Fig. 8 Evolution of Bl/B^ for the LBE run, where Bl is the longitudinal 
component of the magnetic field. The higher initial value is due to the 
unequal physical lengths in the x- and y-directions, associated to the use of 
the hexagonal lattice. 
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